function R = varpe_features(image)
R = [processOne(image, 1), ...
    processOne(image, 2), ...
    processOne(image, 3), ...
    processOne(image, 4), ...
    processOne(image, 5)];
end
   
function result = lbpvar(image, radius, neighbors)
image=double(image);

spoints=zeros(neighbors,2);

% Angle step.
a = 2*pi/neighbors;
for i = 1:neighbors
    spoints(i,1) = -radius*sin((i-1)*a);
    spoints(i,2) = radius*cos((i-1)*a);
end

% Determine the dimensions of the input image.
[ysize xsize] = size(image);

miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));

% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;

% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));

% Minimum allowed size for the input image depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
    error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end

% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;

% Fill the center pixel matrix C.
C = image(origy:origy+dy,origx:origx+dx);
d_C = double(C);

bins = 2^neighbors;

% Initialize the result matrix with zeros.
result=zeros(dy+1,dx+1);

%Compute the LBP code image
S = zeros(dy+1, dx+1);
for i = 1:neighbors
    y = spoints(i,1)+origy;
    x = spoints(i,2)+origx;
    % Calculate floors, ceils and rounds for the x and y.
    fy = floor(y); cy = ceil(y); ry = round(y);
    fx = floor(x); cx = ceil(x); rx = round(x);
    % Check if interpolation is needed.
    if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
        % Interpolation is not needed, use original datatypes
        N = image(ry:ry+dy,rx:rx+dx);
    else
        % Interpolation needed, use double type images
        ty = y - fy;
        tx = x - fx;
        
        % Calculate the interpolation weights.
        w1 = (1 - tx) * (1 - ty);
        w2 =      tx  * (1 - ty);
        w3 = (1 - tx) *      ty ;
        w4 =      tx  *      ty ;
        % Compute interpolated pixel values
        N = w1*image(fy:fy+dy,fx:fx+dx) + w2*image(fy:fy+dy,cx:cx+dx) + ...
            w3*image(cy:cy+dy,fx:fx+dx) + w4*image(cy:cy+dy,cx:cx+dx);
    end
    
    S = S + N;
end
% calculate mean value
S = S / neighbors;

for i = 1:neighbors
    y = spoints(i,1)+origy;
    x = spoints(i,2)+origx;
    % Calculate floors, ceils and rounds for the x and y.
    fy = floor(y); cy = ceil(y); ry = round(y);
    fx = floor(x); cx = ceil(x); rx = round(x);
    % Check if interpolation is needed.
    if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
        % Interpolation is not needed, use original datatypes
        N = image(ry:ry+dy,rx:rx+dx);
    else
        % Interpolation needed, use double type images
        ty = y - fy;
        tx = x - fx;
        
        % Calculate the interpolation weights.
        w1 = (1 - tx) * (1 - ty);
        w2 =      tx  * (1 - ty);
        w3 = (1 - tx) *      ty ;
        w4 =      tx  *      ty ;
        % Compute interpolated pixel values
        N = w1*image(fy:fy+dy,fx:fx+dx) + w2*image(fy:fy+dy,cx:cx+dx) + ...
            w3*image(cy:cy+dy,fx:fx+dx) + w4*image(cy:cy+dy,cx:cx+dx);
    end
    result = result + abs(N-S) / neighbors;
end
%MAPPING = getmapping(8, 'u2');
result = lpq(result, 2*radius+1, 0, 1, 'nh');
end

function f = processOne(image, k)
r = compute_residual_image_succ(image, k);
f = [lbpvar(r, 1, 8), ...
    lbpvar(r, 2, 8), ...
    lbpvar(r, 3, 8)];
end
